## Deserving Government Assistance? 
## Public Support for Aid to Struggling Firms and Workers
#############################################################
## Chris Witko
## Tobias Heinrich
#############################################################
## 3) Estimate average marginal component effects
#############################################################


## Load prepped survey data
###########################
data <- read.csv(file="output/Prepped_survey_data.csv")
CCES <- read.csv(file="output/Prepped_CCES_data.csv")
store_results <- vector("list", 2)


## AMCEs; Firms
###############
this_data <- subset(data, Which_experiment == "Firm")
form <- Choice ~ 1 + I(F_industry == 'Arts, entertainment, recreation') +
  + I(F_industry == 'Manufacturing') + I(F_industry == 'Finance and insurance') +
  + I(F_industry == 'Health care') + I(F_industry == 'Real estate and rental') +
  + I(F_industry == 'Retail') + I(F_industry == 'Transportation and warehousing') +
  + I(F_industry == 'Travel and food services') + I(F_industry == 'Wholesale trade') +
  + I(F_typicalworkerage == '20s') + I(F_typicalworkerage == '30s') + I(F_typicalworkerage == '50s') +
  + I(F_typicalworkerage == '60s') +
  + I(F_numberofworkers == 'Less than 5') + I(F_numberofworkers == 'Between 10 and 19') +
  + I(F_numberofworkers == 'Between 20 and 99') + I(F_numberofworkers == 'Between 100 and 499') +
  + I(F_numberofworkers == '500 or more') +
  + I(F_typicalworkersalary == '$20,000') + I(F_typicalworkersalary == '$160,000') +
  + I(F_typicalworkereducation == '4-year university') +
  + I(F_ceopay == '$10 million') + I(F_ceopay == '$200,000') +
  + I(F_typicalworkerraceethnicity == 'Black') +
  + I(F_typicalworkerraceethnicity == 'Hispanic') + I(F_typicalworkerraceethnicity == 'Asian') +
  + I(F_typicalworkerraceethnicity == 'Mixed') + I(F_typicalworkerraceethnicity == 'Other') +
  + I(F_typicalworkerraceethnicity == 'Native American') +
  + I(F_reasonforstruggle == 'Competition due to superior use of technology (like robots, artificial intelligence)') +
  + I(F_reasonforstruggle == 'Foreign competition due to lower labor cost') +
  + I(F_reasonforstruggle == 'Mismanagement') +
  + I(F_reasonforstruggle == 'Operations closed due to pandemic public health restrictions') 

## Entropy balancing
CCES$id <- sample(1:10, size=nrow(CCES), replace=TRUE)
eb_form <- id ~ 1 + Gender_female + Age + Education_uni + I(Party == "Democrat") + I(Party == "Republican") +
  + I(Job == "Full-time") + I(Job == "Unemployed")
X_cces <- model.matrix(eb_form, CCES)[, -1]
X_data <- model.matrix(eb_form, this_data)[, -1]
X_eb <- rbind(X_data, X_cces)

eb <- ebalance(Treatment=c(rep(0, nrow(this_data)), rep(1, nrow(X_cces))), X=X_eb) 
eb <- ebalance.trim(eb)
w <- eb$w / sum(eb$w) * nrow(this_data)


## Estimate models, weighted/ unweighted  
coefs <- c()
all_amce <- c()
for(j in 1:2)
{
  if(j == 1) mod <- lm(form, data=this_data, weights=w)
  if(j == 2) mod <- lm(form, data=this_data)
  beta <- coef(mod)
  VCV <- vcovBS(x=mod, cluster=this_data$id, R=250)
  store_coef <- rmvnorm(n=1000, mean=beta, sigma=VCV)
  store_sigma <- (summary(mod)$sigma)
  
  ## Summarize for table
  coefs <- rbind(coefs, ddply(.data=melt(cbind(store_coef, store_sigma)), .variables="Var2",
                              .fun=function(x) data.frame(Weights=ifelse(j == 1, "Yes", "No"), What="Firms",
                                                          Mean=mean(x$value), SD=sd(x$value),
                                                          Q025=quantile(x$value, 0.025),
                                                          Q975=quantile(x$value, 0.975))))
  
  ## Put together AMCEs
  zeros <- matrix(0, nrow(store_coef), 16)
  colnames(zeros) <- c(paste0(c("I(F_industry == \"Construction\")",
                                "I(F_typicalworkerage == \"40s\")",
                                "I(F_numberofworkers == \"Between 5 and 9\")",
                                "I(F_typicalworkersalary == \"$70,000\")",
                                "I(F_typicalworkereducation == \"Less than 4-year university\")",
                                "I(F_ceopay == \"$2 million\")", "I(F_typicalworkerraceethnicity == \"White\")",
                                "I(F_reasonforstruggle == \"Economic recession\")"), "TRUE"), 
                       "Industry", "Typical worker age", "Number of workers", "Typical salary", "Typical education", "CEO pay", "Typical race/ ethnicity",
                       "Reason for struggle")
  tmp <- cbind(store_coef[, -1], zeros)
  
  tmp <- adply(.data=tmp, .margins=2, 
               .fun=function(x) data.frame(Weights=ifelse(j == 1, "Yes", "No"), What="Firms", Mean=mean(x), SD=sd(x), PPP=mean(x > 0),
                                           Q025=quantile(x, 0.025), Q975=quantile(x, 0.975)))
  tmp$Header <- ifelse(tmp$X1 %in% c("Industry", "Typical worker age", "Number of workers",
                                     "Typical salary", "Typical education", "CEO pay", "Typical race/ ethnicity",
                                     "Reason for struggle"), 1, 0)
  all_amce <- rbind(all_amce, tmp)
}

## Storage
store_results[[1]] <- list(AMCE=all_amce, Coefs=coefs)


## AMCEs; Workers
#################
this_data <- subset(data, Which_experiment == "Worker")
form <- Choice ~ 1 + I(W_salary == "$20,000") + I(W_salary == "$160,000") +
  + I(W_raceethnicity == 'Black') +
  + I(W_raceethnicity == 'Hispanic') + I(W_raceethnicity == 'Asian') +
  + I(W_raceethnicity == 'Mixed') + I(W_raceethnicity == 'Other') +
  + I(W_raceethnicity == 'Native American') +
  + I(W_education == '4-year university') + I(W_gender == "Female") +
  + I(W_age == '20s') + I(W_age == '30s') + I(W_age == '50s') + I(W_age == '60s') +
  + I(W_reasonforjobloss == "Company was closed due to pandemic public health restrictions") +
  + I(W_reasonforjobloss == "Fired for poor performance") +
  + I(W_reasonforjobloss == "Job was replaced by technology (like robot, artificial intelligence)") +
  + I(W_reasonforjobloss == "Laid off due to company mismanagement") +
  + I(W_reasonforjobloss == "Laid off due to domestic competition") +
  + I(W_reasonforjobloss == "Laid off due to foreign competition") +
  + I(W_reasonforjobloss == "Operations were moved to a foreign country")


## Entropy balancing
CCES$id <- sample(1:10, size=nrow(CCES), replace=TRUE)
X_cces <- model.matrix(eb_form, CCES)[, -1]
X_data <- model.matrix(eb_form, this_data)[, -1]
X_eb <- rbind(X_data, X_cces)
eb <- ebalance(Treatment=c(rep(0, nrow(this_data)), rep(1, nrow(X_cces))), X=X_eb) 
eb <- ebalance.trim(eb)
w <- eb$w / sum(eb$w) * nrow(this_data)
this_data$w <- w

## Estimate models, weighted/ unweighted  
coefs <- c()
all_amce <- c()
for(j in 1:2)
{
  if(j == 1) mod <- lm(form, data=this_data, weights=w)
  if(j == 2) mod <- lm(form, data=this_data)
  beta <- coef(mod)
  VCV <- vcovBS(x=mod, cluster=this_data$id, R=250)
  store_coef <- rmvnorm(n=1000, mean=beta, sigma=VCV)
  store_sigma <- summary(mod)$sigma
  
  ## Summarize for table
  coefs <- rbind(coefs, ddply(.data=melt(cbind(store_coef, store_sigma)), .variables="Var2",
                              .fun=function(x) data.frame(Weights=ifelse(j == 1, "Yes", "No"), What="Workers", Mean=mean(x$value), SD=sd(x$value),
                                                          Q025=quantile(x$value, 0.025),
                                                          Q975=quantile(x$value, 0.975))))
  
  
  ## Calculate AMCEs
  zeros <- matrix(0, nrow(store_coef), 12)
  colnames(zeros) <- c(paste0(c("I(W_age == \"40s\")", "I(W_gender == \"Male\")", 
                                "I(W_salary == \"$70,000\")",
                                "I(W_education == \"Less than 4-year university\")",
                                "I(W_raceethnicity == \"White\")",
                                "I(W_reasonforjobloss == \"Laid off due to economic recession\")"), "TRUE"),
                       "Age", "Gender", "Salary", "Education", "Race/ ethnicity",
                       "Reason for job loss")
  tmp <- cbind(store_coef[, -1], zeros)
  tmp <- adply(.data=tmp, .margins=2, 
               .fun=function(x) data.frame(Weights=ifelse(j == 1, "Yes", "No"), What="Workers", Mean=mean(x), SD=sd(x), PPP=mean(x > 0),
                                           Q025=quantile(x, 0.025), Q975=quantile(x, 0.975)))
  tmp$Header <- ifelse(tmp$X1 %in% c("Age", "Gender", "Salary", "Education", "Race/ ethnicity",
                                     "Reason for job loss"), 1, 0)
  
  all_amce <- rbind(all_amce, tmp)
}

## Storage
store_results[[2]] <- list(AMCE=all_amce, Coefs=coefs)


## Save
#######
all_coefs <- rbind(store_results[[1]]$Coefs, store_results[[2]]$Coefs)
all_amces <- rbind(store_results[[1]]$AMCE, store_results[[2]]$AMCE)
write.csv(x=all_amces, file="output/All-AMCEs.csv")
write.csv(x=all_coefs, file="output/All-Coefs.csv")


## Rope-ladder plots of AMCEs
#############################
all_amces$Sig <- ifelse(all_amces$PPP > 0.975 | all_amces$PPP <= 0.025, "Sig", "Insig")
all_amces$Sig[all_amces$Header == 1] <- "Other"
all_amces$Sig[all_amces$Q975 == 0] <- "Other"
for(j in 1:2)
{
  for(i in 1:2)
  {
    if(i == 1)
    {
      tmp <- subset(all_amces, What == "Workers" & Weights == ifelse(j == 1, "Yes", "No"))
      tmp$PrettyLabel <- mapvalues(x=tmp$X1, 
                                   from=c(paste0("I(W_raceethnicity == \"", c("White", "Black", "Hispanic", "Asian", "Native American", 
                                                                              "Mixed", "Other"), "\")TRUE"),
                                          paste0("I(W_salary == \"", c("$20,000", "$70,000", "$160,000"), "\")TRUE"),
                                          paste0("I(W_education == \"", c("Less than 4-year university", "4-year university"), "\")TRUE"),
                                          paste0("I(W_gender == \"", c("Male", "Female"), "\")TRUE"),
                                          paste0("I(W_age == \"", c("20s", "30s", "40s", "50s", "60s"), "\")TRUE"),
                                          paste0("I(W_reasonforjobloss == \"", c("Laid off due to economic recession",
                                                                                 "Fired for poor performance",
                                                                                 "Company was closed due to pandemic public health restrictions",
                                                                                 "Job was replaced by technology (like robot, artificial intelligence)",
                                                                                 "Laid off due to company mismanagement",
                                                                                 "Laid off due to domestic competition",
                                                                                 "Laid off due to foreign competition",
                                                                                 "Operations were moved to a foreign country"), "\")TRUE"),
                                          "Age", "Gender", "Salary", "Education", "Race/ ethnicity", "Reason for job loss"),
                                   to=c(paste0("       ", c("White", "Black", "Hispanic", "Asian", "Native American",
                                                            "Mixed", "Other",
                                                            "$20,000", "$70,000", "$160,000",
                                                            "Less than 4-year university", "4-year university",
                                                            "Male", "Female",
                                                            "20s", "30s", "40s", "50s", "60s",
                                                            "Economic recession", "Poor performance", "Pandemic public health", "Technology",
                                                            "Company mismanagement",  "Domestic competition", "Foreign competition", "Offshoring")),
                                        "Age", "Gender", "Salary", "Education", "Race/ ethnicity", "Reason for job loss"))
      tmp$PrettyLabel <- factor(tmp$PrettyLabel, 
                                levels=rev(c("Salary", paste0("       ", c("$20,000", "$70,000", "$160,000")),
                                             "Reason for job loss", paste0("       ", c("Economic recession", "Poor performance", "Pandemic public health", "Technology",
                                                                                        "Company mismanagement",  "Domestic competition", "Foreign competition", "Offshoring")),
                                             "Race/ ethnicity", paste0("       ", c("White", "Black", "Hispanic", "Asian", "Native American",
                                                                                    "Mixed", "Other")),
                                             "Age", paste0("       ", c("20s", "30s", "40s", "50s", "60s")),
                                             "Education", paste0("       ", c("Less than 4-year university", "4-year university")),
                                             "Gender", paste0("       ", c("Female", "Male")))))
    }
    if(i == 2)
    {
      tmp <- subset(all_amces, What == "Firms" & Weights == ifelse(j == 1, "Yes", "No"))
      tmp$PrettyLabel <- mapvalues(x=tmp$X1, 
                                   from=c(paste0("I(F_industry == \"", c("Arts, entertainment, recreation", "Manufacturing",
                                                                         "Finance and insurance", "Health care",
                                                                         "Real estate and rental", "Retail",
                                                                         "Transportation and warehousing",
                                                                         "Travel and food services", "Wholesale trade",
                                                                         "Construction"), "\")TRUE"),
                                          paste0("I(F_typicalworkerage == \"", c("20s", "30s", "40s", "50s", "60s"), "\")TRUE"),                                        
                                          paste0("I(F_numberofworkers == \"", c("Less than 5", "Between 5 and 9",
                                                                                "Between 10 and 19", 
                                                                                "Between 20 and 99",
                                                                                "Between 100 and 499",
                                                                                "500 or more"), "\")TRUE"),
                                          paste0("I(F_typicalworkereducation == \"", c("Less than 4-year university", "4-year university"), "\")TRUE"),
                                          paste0("I(F_typicalworkersalary == \"", c("$20,000", "$70,000", "$160,000"), "\")TRUE"),
                                          paste0("I(F_ceopay == \"", c("$200,000", "$2 million", "$10 million"), "\")TRUE"),
                                          paste0("I(F_typicalworkerraceethnicity == \"", c("White", "Black", "Hispanic", "Asian", "Native American",
                                                                                           "Mixed", "Other"), "\")TRUE"),
                                          paste0("I(F_reasonforstruggle == \"", c("Competition due to superior use of technology (like robots, artificial intelligence)",
                                                                                  "Foreign competition due to lower labor cost",
                                                                                  "Mismanagement", "Operations closed due to pandemic public health restrictions",
                                                                                  "Economic recession"), "\")TRUE"),
                                          "Industry", "Typical worker age", "Number of workers", "Typical salary",
                                          "Typical education", "CEO pay", "Typical race/ ethnicity", "Reason for struggle"),
                                   to=c(paste0("       ", c("Arts, entertainment, recreation", "Manufacturing",
                                                            "Finance and insurance", "Health care", "Real estate and rental", "Retail",
                                                            "Transportation and warehousing", "Travel and food services", "Wholesale trade",
                                                            "Construction", 
                                                            "20s", "30s", "40s", "50s", "60s",
                                                            "Less than 5", "Between 5 and 9",
                                                            "Between 10 and 19",  "Between 20 and 99", "Between 100 and 499",
                                                            "500 or more", "Less than 4-year university", "4-year university", 
                                                            "$20,000", "$70,000", "$160,000", "$200,000", "$2 million", "$10 million",
                                                            "White", "Black", "Hispanic", "Asian", "Native American",
                                                            "Mixed", "Other",
                                                            "Technology", "Foreign competition", "Mismanagement", "Pandemic public health",
                                                            "Economic recession")),
                                        "Industry", "Typical age", "Number of workers", "Typical salary",
                                        "Typical education", "CEO pay", "Typical race/ ethnicity", "Reason for struggle"))
      
      tmp$PrettyLabel <- factor(tmp$PrettyLabel, 
                                levels=rev(c("Typical salary", 
                                             paste0("       ", c("$20,000", "$70,000", "$160,000")),
                                             "CEO pay", 
                                             paste0("       ", c("$200,000", "$2 million", "$10 million")),
                                             "Reason for struggle", paste0("       ", c("Pandemic public health", "Economic recession", "Foreign competition",
                                                                                        "Technology", "Foreign competitiont", "Mismanagement")),
                                             "Industry", 
                                             paste0("       ", c("Health care", "Construction", "Manufacturing", "Travel and food services", 
                                                                 "Transportation and warehousing", "Retail",
                                                                 "Wholesale trade", "Arts, entertainment, recreation", 
                                                                 "Finance and insurance", "Real estate and rental")),
                                             "Typical race/ ethnicity", 
                                             paste0("       ", c("White", "Black", "Hispanic", "Asian", "Native American",
                                                                 "Mixed", "Other")), 
                                             "Typical age", paste0("       ", c("20s", "30s", "40s", "50s", "60s")),
                                             "Number of workers", 
                                             paste0("       ", c("Less than 5", "Between 5 and 9",
                                                                 "Between 10 and 19",  "Between 20 and 99", "Between 100 and 499",
                                                                 "500 or more")), 
                                             "Typical education", 
                                             paste0("       ", c("Less than 4-year university", "4-year university"))
                                )))
    }
    
    ## Rescale
    tmp[, c("Mean", "Q025", "Q975")] <- 100 * tmp[, c("Mean", "Q025", "Q975")]
    tmp$Mean_label <- format(round(tmp$Mean, di=0), nsmall=0)
    tmp$Base <- 0
    tmp$Base[tmp$Header == 0 & tmp$Q025 == 0] <- 1
    
    for(k in 1:2)
    {
      tmp2 <- tmp
      if(k == 1)
      {
        tmp2$Relevant <- 0
        tmp2$Relevant[str_detect(string=tmp2$X1, pattern="W_reasonforjobloss")] <- 1
        tmp2$Relevant[str_detect(string=tmp2$X1, pattern="W_salary")] <- 1
        tmp2$Relevant[str_detect(string=tmp2$X1, pattern="Salary")] <- 1
        tmp2$Relevant[str_detect(string=tmp2$X1, pattern="Reason for job loss")] <- 1
        tmp2$Relevant[str_detect(string=tmp2$X1, pattern="F_ceopay")] <- 1
        tmp2$Relevant[str_detect(string=tmp2$X1, pattern="F_typicalworkersalary")] <- 1
        tmp2$Relevant[str_detect(string=tmp2$X1, pattern="F_reasonforstruggle")] <- 1
        tmp2$Relevant[str_detect(string=tmp2$X1, pattern="CEO pay")] <- 1
        tmp2$Relevant[str_detect(string=tmp2$X1, pattern="Typical salary")] <- 1
        tmp2$Relevant[str_detect(string=tmp2$X1, pattern="Reason for struggle")] <- 1
        tmp2 <- subset(tmp2, Relevant == 1)
      }

      g <- ggplot(data=tmp2, aes(x=PrettyLabel, y=Mean, colour=factor(Header)))
      g <- g + geom_hline(yintercept=0, linewidth=0.8)
      g <- g + geom_point(size=0.00002, alpha=1/1000)
      g <- g + geom_point(data=subset(tmp2, Base == 1), size=3.5, colour="black")
      g <- g + scale_colour_manual(values=c("black", "white"))
      g <- g + coord_flip() + theme_minimal()
      g <- g + geom_linerange(data=subset(tmp2, Sig == "Sig"),
                              aes(x=PrettyLabel, ymin=Q025, ymax=Q975), colour="black", linewidth=.99)
      g <- g + geom_linerange(data=subset(tmp2, Sig == "Insig"),
                              aes(x=PrettyLabel, ymin=Q025, ymax=Q975), colour="gray60", linewidth=.99)
      g <- g + geom_point(data=subset(tmp2, Sig == "Insig"), size=3, colour="gray60")
      g <- g + geom_label(data=subset(tmp2, Sig == "Sig"), aes(label=Mean_label), size=3.6)
      g <- g + theme(axis.text.y=element_text(hjust=0, size=rel(1.34)),
                     axis.text.x=element_text(size=rel(1.34)),
                     legend.position="none") 
      g <- g + xlab("") + ylab("Average marginal component effect")
      if(i == 1) ggsave(paste0("output/AMCE-Workers-", ifelse(j == 1, "weighted", "unweighted"), "-", ifelse(k == 1, "short", "full"),
                               ".pdf"), width=11, height=9.3 + ifelse(k == 1, -3, 0),
                        device=cairo_pdf)
      if(i == 2) ggsave(paste0("output/AMCE-Firms-", ifelse(j == 1, "weighted", "unweighted"), "-", ifelse(k == 1, "short", "full"),
                               ".pdf"), width=11, height=9.3 + ifelse(k == 1, -3, 0),
                        device=cairo_pdf)
    }
    rm(tmp)
  }
}


## Garbage collection
#####################
rm(list=ls())



